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Abstract The fractional birth and the fractional death 
processes are more desirable in practice than their clas- 
sical counterparts as they naturally provide greater flex- 
ibility in modeling growing and decreasing systems. In 
this paper, we propose formal parameter estimation 
procedures for the fractional Yule, the fractional lin- 
ear death, and the fractional sublinear death processes. 
The methods use all available data possible, are com- 
putationally simple and asymptotically unbiased. The 
procedures exploited the natural structure of the ran- 
dom inter-birth and inter-death times that are known to 
be independent but are not identically distributed. We 
also showed how these methods can be applied to cer- 
tain models with more general birth and death rates. 
The computational tests showed favorable results for 
our proposed methods even with relatively small sample 
sizes. The proposed methods are also illustrated using 
the branching times of the plethodontid salamanders 
data of Highton and Larson (1979). 
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1 Introduction 

Recently, generalizations of the classical birth and death 
processes have been developed using the techniques of 
fractional calculus. These are called the fractional birth 
(Uchaikin et al. 2008; Orsingher and Polito 2010; Ca- 
hoy and Polito 2012) and the fractional death (Ors- 
ingher et al. 2010) processes, correspondingly. A ma- 
jor advantage of these models over their classical coun- 
terparts is that they can capture both Markovian and 
non-Markovian structures of a growing or decreasing 
system. 

When the birth and death rates are both linear, they 
are then called the fractional linear birth or fractional 
Yule or Yule-Furry process (fYp) and fractional linear 
death process, respectively. The classical linear birth 
or Yule process has been widely used to model various 
stochastic systems such as cosmic showers in physics 
and epidemics in biology to name a few (see e.g., Nee 
et al. 1994a; Aldous 2001; Nee 2001; Paradis 2012). 
Note also that the fractional linear birth process was 
partially investigated by Uchaikin et al. (2008) using 
the Ricmann-Liouville derivative operator but was con- 
tinued and generalized by Orsingher and Polito (2010) 
using the Caputo derivative. The inter-birth time dis- 
tribution, which provided a way to simulate the fYp 
was derived in Cahoy and Polito (2012). With this, 
we adopt the fYp from Orsingher and Polito (2010). 
In addition, the definition of the fractional linear and 
fractional sublinear death processes are taken from Ors- 
ingher et al. (2010). 

For completeness, we first enumerate some proper- 
ties of the fractional Yule (with one progenitor) and 
the fractional linear death (with initial population size 
no > 1) processes, which will be used in the subse- 
quent discussions. Table 1.1 below shows the probabil- 
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ity P(t) of no event (no birth or no death) at time t, 
the state probability mass function Pi(t) or the proba- 
bility of having (i — 1) births or (n Q — i) deaths by time 
t, the probability density function fi(t) of the indepen- 
dent but non-identically distributed random inter-event 
times, the mean, and the variance of the fractional Yule 
and the fractional linear death processes. Note that the 
fractional Yule and the fractional linear death processes 
have the parameters A > and fi > as the birth and 
death intensities, correspondingly. 



Note that 



Es,p {%) 



OC 



(1.1) 



is the Mittag-Leffler function. 1 

In this article, we propose regression-based proce- 
dures to estimate the parameters of the fractional lin- 
ear birth, the fractional linear death, and the fractional 
sublinear death processes. The rest of the paper is or- 
ganized as follows. In Section 2, the specific functional 
forms of the inter-death time distributions and the vari- 
ances of the fractional linear and sublinear death pro- 
cesses are obtained. These results allowed us to apply 
our methods to these processes. Section 3 introduces 
the proposed method using the fractional Yule, the frac- 
tional linear death, and the fractional sublinear death 
processes as examples. The section also shows some ex- 
tensions of the procedures to certain models. Section 
4 contains the empirical test results and the real-data 
application of the proposed methods for the case of the 
fYp only as similar inference procedures can be applied 
to the fractional linear and fractional sublinear death 
processes. The summary and extensions of our study 
are given in Section 5. 



2 More properties of the fractional linear and 
fractional sublinear death processes 

We now derive some properties which will permit us to 
apply the proposed estimation procedures to the frac- 
tional linear death and the fractional sublinear death 
processes. More specifically, the theorems below showed 
that the inter-death times for both the fractional linear 
and sublinear death processes are Mittag-Leffler dis- 
tributed. The variances of both processes are also de- 
rived. 

1 Note: The entries with (**) are new results and are de- 
rived in Section 2. 



Theorem 2.1 The inter- death time T£ of the fractional 
linear death process {M u (t),t > 0} with death rate in- 
tensity \i > 0, and n Q G N initial individuals are inde- 
pendent but are non-identically distributed with proba- 
bility density function 

Pr{2£ € dt}/dt = n(no - fc)f- 1 S I/)V (-^(n - k)t v ), 

where k = 0, 1, . . . , no — 1, and T% is the random time 
separating the kth and (k + l)th death. 

Proof We prove the theorem by induction. When k = 
we obtain 



(2.1) 



Pr{T " <t}=Pr{M"(t) < n } 

= l-Pr{M v (t)=n } 

= 1 - P(t) (see Table 1.1) 

= 1 - E v ,i{-^n Q t v ). 

Therefore 

Pv{T% G dt}/dt = j Px{T% < t} (2.2) 

= fj,not u ~ Ev^-^nat"). 
For k = 1 we observe 

Pr{T l/ +T? G dt}/dt (2.3) 
= jVv{T^ + T^<t] 

= ±p v {M v {t) <n -l} 
at 

= j t [l- Pr{M v (t) = n } - Pr{M"(t) = n - 1}] . 



Using Table 1.1 we get 

Pr{T w + T" G dt}/dt 

= - ^iE v x(-imot v ) 
dt 



(2.4) 



- "77 [noE v i(-(riQ - 1)/J,t v ) - n E v ,x(-noiit v )} 
dt 

= pLntf v ~ 1 E v . y (-iin§t v ) 

+ n Q (n - l)/xt I/_1 ^, r/ (-ju(no - l)t") 

- n\iit v ^ 1 E vy (~[inQt v ) 

= n (n - l)^^ 1 [E Vt „(-(j,(n - l)t") 
-E UiU (-fj,n t v )] . 

To check the preceding results, we can obtain the Laplace 
transform as 



/>oo 

/ e~ wt Pr{2£ + Ty G dt} 
Jo 

n (n Q - l)fi _ n (n - l)/z 
w v + fx(no — 1) w v + /jhq 



(2.5) 



Parameter estimation for fractional birth and fractional death processes 



3 



Table 1.1 Known properties of fractional Yule (JV" (t)) and linear death (M"(t)) processes. 





fractional Yule process 


fractional linear death process 


Pit) 
Pi(t) 

Mean 
Variance 


£„,l(-At") 

\it v - 1 E v ^{-\it v ), i>l 
E v>1 (At") 

2E„ A (2XF) - E uA (At") - (E„ A (At")) 2 


(?) e;^o 1 {^(-lyE^i-ii+j^n, o<i< no 

H(n - i)t"- 1 E v>v (-p,{na - i)t"), < i < n - 1 (**) 
n E v ,i (-fit") 

n (n - l)E uA (-2^) + n E Uil (-fit") - n% (E U:1 (-^t v )) 2 (**) 



fj,n /Lt(n - 1) 



w v + fj.no w v + jLt(no — 1) 

/■OO /"OO 

= / e- ws Pt{T£ G ds} / e~ wy Pt{T? G ds} 
Jo Jo 

/•OO f-t 

= / e~ wt / Pr{2t G d(t - s)} Pr{T» G ds} 
Jo Jo 

/"OO />oo 

= ^ Pr{2£ € ds} J e- zt Pv{T^ G d(t - s)}, 

which is just a convolution of two independent variables 
Tq and T" . For a general k it is sufficient to note that 

Pr{T£ + ...+Tj;edt} (2.6) 

= / Pr{T£ G d(t - a)} Pr{TJ + • • • + T£_ x G ds}. 
Jo 

By exploiting again the Laplace transform and writing 
®k = T a +---+ T k> we have 

pOO 

/ e~ wt Pr{££ G dt} (2.7) 
Jo 

t-OO f-t 

= / e~ wt / Pr{T£ G d(t - «)} Pr{D^_ 1 G <fe} 
Jo Jo 

poo poo 

= J PriS)^ G ds} J e -* Pr{T£ G d(t - s)} 

pOO pOO 

- / e~ ws Pr{©j;-i G ^} / e""" Pr{T£ G dy} 
Jo •/ 

TT / e- ws Pr{T/ G ds} 

,,-nJO 



J=0 
fc 

n 

J=0 



n ( n o - J) 
w v + /i(n - j) ' 



We now determine the variance of the fractional lin- 
ear death process {M"(t),t > 0}. Consider equation 
(1.6) of Orsinghcr et al. (2010). That is, 

f £,p%(t) = fl(k + l)pl +1 {t) - nhpl{t), < fc < n , 



^(0) 



1, k = n , 
0, < k < no. 



It is then straightforward to arrive at 

(c^ G » M = -n{u- l)&G"(t*,t), 
\G v (u,0) = u n °, 



(2-9) 



where G v (u,t) — J2k=o uk Pk(t) * s the probability gen- 
erating function of the fractional linear death process. 
This in turn leads to 



• r H(t) = -2uH{t), 



H(0) = n (n - 1), 



(2.10) 



where H(t) = E(M V '(t)(M u '(t) - 1)) is the second fac- 
torial moment. The solution to (2.10) reads 

H(t) = n (n - l)£? J/)1 (-2/if) ) (2.11) 

and the variance can be immediately obtained as 

VarM v (t) = H(t) + EM" (t) - (EM"(t) f (2.12) 
= no{no-l)E V;l {-2nt v ) 
+ noE^i-ut") - n 2 Q (E uA (-ut»)) 2 . 

Note that the above expression, when v = 1, simplifies 
to the variance of the classical linear death process, i.e. 



VarM 1 ^) = n e-' 1 *(l - e~^). 



(2.13) 



(2.8) 



Below is the algorithm to generate a typical sam- 
ple path of a fractional linear death process in Figure 
2.1. Note that there are several sub- algorithms to gen- 
erate the inter-death times T!j"s that are available in 
the literature (see e.g., Cahoy and Polito 2012). 

Algorithm: 

Step 1. Let k — and the population size equal no- 
Step 2. Simulate T% , and let the A;th death time be 

®l = T» + Tf + T$ + ■ ■ ■ + T». 
Step 3. Set the population size no — k, and k = k + 1. 
Step 4. Repeat Steps 2-3 for fc = 1, . . . , n — 1. 
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of Orsingher et al. (2010). Recall that 
d 2 



du 2 



= E[2T(t) (£W(t)-l)] 
= H{t). 



(2.14) 



Then 

d v 



H(t) = -2/i(n + l) (E*m u (t) + Pr{art"(i) = 0} - 1) 




Fig. 2.1 Sample paths of the classical linear death process 
(top) and the fractional linear death process (bottom) in the 
interval with parameters (v, A) = (0.75, 1) and initial popula- 
tion size no = 40. 



It can be gleaned from Figure 2.1 that the sam- 
ple path of the fractional linear death process (bot- 
tom) seems to decay faster at small times but is slower 
for large times than its classical counterpart. The fig- 
ure also indicates that it is capable of producing death 
bursts especially at early stages (corresponding to small 
times). 

The inter-death time distribution for the fractional 
sublinear death process can be easily deduced (whose 
proof follows from the previous result and is omitted) 
from the preceding theorem as follows. 

Theorem 2.2 The fractional sublinear death process 
{fSt v (t), t > 0}, with death intensity rate \i > 0, and 
no € N initial individuals has the following probability 
density function of the inter-death times 's 

Pr{1£ € dt}/dt = fj,(k + ly^E^-^k + l)t u ), 



with k = 0, 1, 



,n 



1, where T£ is the random time 



separating the kth and (k + l)th death. 

The variance of the fractional sublinear death pro- 
cess can be determined by considering equation (3.45) 



2/iif(t) 



(2.15) 



\k=l 



2/x(n + l) (£(?) (-1)^,1 



£ (^ + 1 1 )(-l) fe+1 ^, 1 (-/i^)^ + 2/*ff(t) 



= -2/i(n +l) 



E 



Lk = l 

2fj,H(t) 



n Q + l 
k + l 



{-\) k E uA (-knt v ) 



k=l 



2^n + l)Y,[ k n ° 1 )(-l) k E v A-^t 



+ 2/j,H(t). 

Using the initial condition H(0) — no(no — 1) and let- 
ting H(w) be the Laplace transform of H(t), we write 



w u H(w)—w v no(riQ — 1) 



(2.16) 



k=l 



Hence, 
H(w) 



k[i 



2fj,H(w) 



(2.17) 



i ^\ w , . \ ■> / no 



k=l 



k„.,v-\ 



n (n Q - 1) 



{w v + k^){w v - 2/x) 



w u - 2fi 



+ 2n(n Q + 1) 



fe=i v 
n (n -l 



1 



1 



w v + kfi w lJ — 2/i 



(-2A*) 



w u -2[i w v - 2^, 



k=l 
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■("o + i)E 



k=l 



,(-1) - 

k + 1 1 v ' w v + k^i 



w u - 2fi 



(l-no)-(no+l)2 



k = l 



k + lj y ' a)" + k[i 



The second factorial moment can be easily shown as 



H(t) = -(n - 1)E V<1 (2^) 



(2.18) 



>o + l)£ 



fc=i 



u 

fc + 1 



(-l^+^xC-fc/iO- 



Thus, the variance simply follows as 
WarWl u {t) = ff(t) + E93T(t) - [£W(t)] 

"0 

= -(n - 1)^,1(2^) + (n + 1) E 



fc=i 



n 
fc + 1 



x(-i) fc+i £;, il (-fc^ l/ ) 



E(?^)(-i) fc+1 ^.i(-^) 



no 

E 

Lfc=l 



n + 1 
fc + 1 



(-l) k+1 E Utl (-k^) 



Note that the algorithm above could be easily adopted 
to simulate sample trajectories of the fractional sublin- 
ear death process. 



expectation on both sides, it can be easily shown that 
the mean and variance (see details in Cahoy et al. 
2010) of the log-transformed i-th random sojourn time 
T- = ln(Tj) of the fYp are 



and 



In (At) 



1 

3^ 



- 7, 

1 

6 



(3.1) 



(3.2) 



respectively, where 7 ~ 0.5772156649 is the Euler - 
Mascheroni's constant. The first two moments above 
therefore suggest that the following simple linear re- 
gression model can be fitted/formulated: 



T t = a Q + ai In i + e, 
where 

-ln(A) 

a = 7, 

v 



1. 



«i 



(3.3) 



(3.4) 



a^t ) ■ The trick used here 



and £i l = JV (fj, e = 0, 

was to factor out the non-identical means of the log- 
transformed random inter-birth or sojourn times, which 
are linear functions of the logarithm of the known fixed 
i. Thus, this leads to studying the widely used simple 
linear regression model (see Montgomery et al. 2006). 

3.1.1 Point estimation 



3 Parameter estimation 

3.1 Estimation for the fractional Yule or linear birth 
process 

We now illustrate our estimation approach for the fYp 
with birth rate Ai, i > 1. Furthermore, assume that a 
sample trajectory of n births corresponding to n ran- 
dom inter-birth times TVs of the fractional linear birth 
process is observed. That is, n independent but are not 
identically distributed random inter-birth times of the 
fractional linear birth process are given. This also in- 
sinuates that only a single datum is obtained from each 
of the n different Mittag-Leffler distributions. This ob- 
servation and the Mittag-Lefner's seemingly complex 
structure pose a computational challenge on how to 
estimate the model parameters more efficiently espe- 
cially for small population sizes. Recall the structural 
representation of the Mittag-Leffler distributed random 
inter-birth time Tj = E l / V S u (see Cahoy and Polito 
2012), where E = exp(At) is independent of S v which is 
a one-sided a + -stable distributed random variable. Ap- 
plying the logarithmic transformation and taking the 



Inverting the least squares (LS) estimators 



a 1 



ELi (lnj-Inl) 



(3.5) 



and So = — a± ■ In i gives the LS-based point estima- 
tors of v and A as 



-1 
Si 



and 

Xis = exp ((a + 7)/ ai) 



(3.6) 



(3.7) 



respectively, where lni = lnj/n, and T' = Tj/r 



3=1 



Equating o\ or a^, in (3.1) with its unbiased estimator 



we get the residual-based point estimators 

' (see Cahoy et al. 2010) 



(3.8) 



V 3 0*2 A 2 + s) 



(i 
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and 



A res = exp {-v rea (a + 7)) 



(3.9) 



(3.10) 



of the model parameters v and A, correspondingly where 
Si = T i — T[, and T[ = 1io + a\\ni. Note that the 
residual-based estimators exploit the residuals to esti- 
mate v rather than the negative inverse of the LS esti- 
mate of the slope a\. 

3.1.2 Interval estimation 

We now develop interval estimators using the large- 
sample properties of the least squares estimators 60 and 
bi above. The following result shows the joint asymp- 
totic behavior of the proposed point estimators of v and 
A for the fYp. 

Theorem 3.1 Let < v < 1 and A > 0. Then 
% s -v' 



A;. — A 



d 



IV[0,ncr 2 C] 



where " — > " denotes convergence in distribution, 

Ci C12 

v C21 C2 

Ci = u 4 s-\ 

Cia = C21 = \v z ((lni + ln(A)) /*), 

C 2 = {v\f (l/n+ (hif 

+21n(A)Ihl+(ln(A)) 2 ) /a), 

and s — ^2 (l n J ~ m *) • 

i=l 

Proof Recall the large-sample normality of the least 
squares estimators clq anda\, i.e., 



where the covariance matrix X is defined as 
^ (l/n -\-hxi /s^ — ]ni/i 
— lni/s s" 1 



S = n<r? 



V 



Recall the multivariate delta method (Ferguson 1996): 
If y/n0 n - [3) — ► JV[0, E] i/iera 

vMff(A0 - flC9)) 4 N[0, g(/3) T £g(/3)] . 



Hence, using the delta method above, f3 n = (ao,ai) T , 
g(J3 n ) = {vi s ,\i s ) T , g([3) = {v,\) T , and the Jacobian 
matrix 



exp((a + 7)/ai)/ai 

l/aj -exp((a +7)/ai) (a + 7) / a i 



we obtain the final expression of the covariance matrix 
by simply substituting back ag = — ln(X)/v— 7 and a\ — 
-1/v. ■ 

Corollary 3.1 Approximate (1 — a)100% confidence 
intervals for v and A can 6e deduced as 



(3.11) 



P is ± z a/2 a e Vf s Vs \ 
and 

\( s ± z a/2 a e vi s \is(l/n 

1 /2 

+ (hiI 2 + 21n(A is )Inl+(ln(A is )) 2 )/ S ) , (3.12) 

respectively, where z a / 2 is the (1 — a/2)th quantile of 
the standard normal distribution and < a < 1. 

We now propose another interval estimators which 
utilize the residual-based estimate of v, and a bootstrap 
technique. It can be inferred from Cahoy ct al. (2010) 
that a residual-based (1 — a) 100% confidence interval 
for v can be 



'■'res i z a/2 



P 2 es (32 - 20P 2 es -9 r 4 es ) 



40rc 



(3.13) 



where z a / 2 is defined above. A residual-based (1— a) 100% 
interval estimate for A can also be 



A res ± z a j 2 



+ K 



£l (l/n 



lni 2 /s 



1/2 



(3.14) 



Since the small-sample performance of v res and the 
residual-based interval estimator in (3.13) have been 
shown to perform well already (see, e.g., Cahoy et al. 
2010), we apply a non-parametric percentile bootstrap 
technique to A res using the fixed-regressor approach 
to obtain a small-sample interval estimator of A. This 
well-known procedure is slightly modified by first di- 
viding each residual e by yT — hi, where hi is the ith 
leverage or the ith diagonal entry in the hat matrix 
before sampling from the transformed residuals. Note 
that the division of yl — hi is simply for correction 
as the true variance of the residual £ is Yar E; = 
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(Tg(l — hi) (see Montgomery ct al. 2006). Hence, the 
bootstrap counterpart of A res is calculated as X* es — 
exp (— v* es (<Zq + 7)) where 9* es used the bootstrapped 
transformed or weighted residuals. A clear advantage of 
the asymptotic-based procedures over the re-sampling- 
based ones is that they are faster to calculate especially 
for large sample sizes. 



3.2 Estimation for the fractional linear and the 
fractional sublinear death processes 

Assuming that a sample trajectory of no deaths cor- 
responding to no random inter-death times T%'s of a 
fractional death process is observed. Following the pro- 
cedure for the fractional linear birth process in the pre- 
ceding subsection, we can estimate the parameters v 
and \i by regressing In (T£) with ln(no — k). That is, we 
fit the following simple linear regression model: 



In(T£) = & + &irn(nn-fe)+£fe, 



(3.15) 



where k = 0, . . . , no — 1, bo = — ln(/x)/V— 7, b\ is given in 

(3.4) of subsection 3.1, and Ek = N ^fi E = 0, <*j^ T „j 

Following the methodology in the preceding subsection, 
we can straightforwardly obtain the corresponding LS- 
based point estimates of v and fj, from (3.6) and (3.7) 
as 



Vis 



and 



-1 

h 



Jlis = exp (ib + 7) / M , 



(3.16) 



(3.17) 



no — 1 



respectively, where ln(no — k) = ln ^"° o ^ , In (T v ) = 
b = \MT^ - h ■ ln(n - k), and 

3=0 

. E^o 1 ln ( T j) (M"o ~ J) - ln(n - k)) 

h= V ■ ( 3 - 18 ) 

E"=o 1 ( ln ( n o - j) - ln(n - fc) 

Furthermore, the LS-based interval estimates for v and 
/i directly follow from (3.11) and (3.12) of Corollary 3.1 
in subsection 3.1.2, correspondingly. Hence, the approx- 
imate (1 — a) 100% for v and 11 can be explicitly written 
as 



vis ± z a/2 a e v l: 



\ 



'n -l 

^ (ln(n - j) - ln(n - k) 

3=0 



and 

Xls ± Z Q / 2 (TePi s A; s ^l/n 



+ ln(n -k) +2 ln(A ifl )ln(n - k) + (ln(A ;s )) 2 /s 



1/2 



correspondingly. On the other hand, the residual-based 
point and interval estimators of v and [i immediately 
follow from subsection 3.1 as well, where do is replaced 
by b Jn (3.10), e k = ln(T£) - ln^f), and hf(2£) = 
6 + bi ln(n — k). 

A similar approach can be done to obtain estimates 
for the fractional sublinear death process. That is, we 
regress In (T£) with ln(fc + 1), k = 0, 1, . . . , no — 1, or fit 
the model 



ln(X£) =c + Cl ln(fc + l)+£ fe , 



(3.19) 



and follow the procedures used for fractional Yule and 
the fractional linear death processes. In general, we sim- 
ply replace A, lni, lni by /i, ln(no — k) or ln(fc + 1), and 
ln(no — k) or ln(/c + 1), accordingly in the methods of 
subsection 3.1 to obtain the parameter estimators for 
the fractional linear death and the fractional sublinear 
death processes. 

3.3 Some Extensions 

Assume that a fractional birth or death process exists 
with rates 9j, j — 1,2,..., n, where the jth inter-event 
time Xj is Mittag-Leffler distributed with parameter 
9j . Then the mean of X- = In (Xj ) is 



In ( 



7- 



(3.20) 



Based on the above mean formulation, we use the model 



X 3 = d + di ■ q(j) + Ej 



(3.21) 



to estimate more forms of the parameters or rates under 
the two cases below. 

Case 1: When\n(9j) = m(9)+q(j) for some appropriate 
known functions m(9) and q(j) of the parameter 9 and 
j € N, correspondingly. 

In this case, the general form of the regression model 
that could be used for estimation is 



x 3 = -\ 1 + 



m(9) 



1 



■ Si- 



(3.22) 



Clearly, d = —(i+m(6)/v),d\ = —1/v, and q(j) is the 
regressor variable. Using v res or —\jd\ and inverting 

the least squares estimate bo, we can compute m(9) and 
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9 sequentially. Note that the explicitness of 9 depends 
on the form of m. 

Example 1: When 6j is linear, i.e., 9j — 8j then 
ln(6>j) = ln(6>) + ln(j), where m(9) = \n(9) and q(j) = 
ln(j), respectively. Note that this parametrization cor- 
responds to the fractional Yule, the fractional linear 
death, and the fractional sublinear death processes. 

Example 2: If 9 = e e+j then ln(9j) = 9 + j, where 
m(9) = 9 and q(j) = j, correspondingly. This suggests 
that d = -(7 + 9/v). 

Case 2: When hi(9j) = m{9) -q{j) for some appropriate 
known functions m(9) and q(j) of the parameter 9 and 
j G N, correspondingly. 

The general form of the regression model in this case 

is 

^=-7" — + (3-23) 

Apparently, do = —7,^1 = —m,(ff)jv, and q(j) is the 
predictor variable. Using v res and inverting the least 
squares estimate d±, we can calculate miff) and 9 suc- 
cessively. 

Example 1: If 6j = 6° then In(^) = ln(0) ■ j, where 
m(9) = hx(9) and q(j) = j, correspondingly. This indi- 
cates that d\ = — ln(0) /v. 

Example 2: When 9j = e e "i then ln(0j) = 9-j, where 
m(9) = 9 and q(j) = j, respectively. This shows that 
d x = -6/v. 

4 Method testing and application 

4.1 Empirical test 

For the sake of reproducibility, we now test our proce- 
dures using the fYp as a particular example as simi- 
lar approach can be carried out for both the fractional 
linear death and the fractional sublinear death pro- 
cesses. In point estimation testing, we evaluated the 
finite-sample properties (unbiasedness and homogene- 
ity) by computing the average and the median absolute 
deviation (MAD) of the estimates using 1000 simula- 
tions for sample sizes n = 100, 500, and 1000. These 
values are shown in Table 4.1 below. The relative fluc- 
tuation (RF=100% xMAD/mean) of vi s decreases from 
19.23% (corresponds to v = 0.1, n = 100) to as little 
as 4.41% (with v = 0.95 and n = 1000). On the other 
hand, the residual-based P res 's RF ranges from 4.41% 
(with v = 0.95 and n = 1000) Jo 1.89% (corresponds 
to v — 0.95, n = 1000). While A; s 's RF improves from 



33.94% (corresponds to A = 0.5, n = 100) to 30.48% 
(with A = 5 and n — 1000), A res 's RF decays faster 
from 55% (A = 1, n = 100) to 25.29% (A = 5 and 
n = 1000). In general, the relative fluctuations of the 
residual-based estimators tend to decay faster than the 
LS-based estimates. They are also less bias than the 
LS-based estimators especially for n < 100. Nonethe- 
less, both the residual- and LS-based point estimators 
are asymptotically unbiased as expected. 



Table 4.2 below shows the averaged lower and upper 
95% confidence bounds using the formulae in Section 
3. These bounds used 1000 simulation runs for each of 
the sample sizes n = 15, 30, 100, and 500. Note that the 
residual-based interval estimator A* es utilized 500 boot- 
strap samples and of = ir 2 (l/(39^ es ) — 1/6) is used to 
estimate the error variance in our LS-based procedures. 
Observe that some of the interval estimates for sample 
sizes n = 15, and n — 30 are omitted as they are unre- 
liable due to the multiple error warnings that showed 
up during the computation process. Moreover, the con- 
vergence of the coverage probabilities to their true lev- 
els for the LS-based method is made faster by using 
the error variance estimate u 2 e — n 2 (l/(3Vf s ) — 1/6). 
From Table 4.2, it is apparent that the residual-based 
interval estimates of v are narrower and are better cen- 
tered around the true parameter values than the least- 
squares' even when the sample size is as large as 500. In 
addition, our simulations showed that the asymptotic 
or non-bootstrapped residual-based interval estimator 
of A gives more sensible results than the LS-based pro- 
cedure for small samples. Nevertheless, the LS-based 
interval estimates for A are more accurately centered 
than the bootstrapped residual-based estimates espe- 
cially for large samples. 



The corresponding coverage probabilities and the 
widths of the interval estimates above with a confi- 
dence level of 95% are displayed in Table 4.3. When 
the sample size n — 15, the residual-based interval esti- 
mators have minimum coverage of 90.1% and 91.1% for 
v = 0.95 and A = 0.1, respectively. When n — 500, the 
bootstrap interval estimator of A has coverage probabil- 
ities which are closer to the true confidence level than 
the LS-based procedure for large values of A. However, 
the LS-based estimator of A has a better coverage than 
the bootstrapped residual-based interval estimator for 
small A values. The residual-based interval estimator for 
A seemed to have slower convergence than the LS-based 
method. Furthermore, the residual-based estimator of 
v outperformed the LS-based method as its coverage 
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Table 4.1 Mean point estimates of and dispersions from the true parameters v and X. 



(i/, A) 


Estimator 


n = 100 
Mean MAD 


n = 500 
Mean MAD 


n = 1000 
Mean MAD 


(0.1, 11 


Vis 
IS res 


0.104 0.020 
0.103 0.008 


0.101 0.008 
0.100 0.004 


0.100 0.006 
0.100 0.003 


Ms 
A res 


3.190 0.665 
1.318 0.725 


1.151 0.407 
1.091 0.408 


1.077 0.318 
1.051 0.322 


(0.25,0.1) 


"is 
l/ res 


0.261 0.048 
0.252 0.022 


0.252 0.021 
0.251 0.010 


0.251 0.014 
0.250 0.007 


A is 
A res 


0.119 0.031 
0.131 0.071 


0.106 0.025 
0.109 0.044 


0.103 0.021 
0.106 0.033 


(0.5,0.5) 


V~l s 

i/ res 


0.521 0.100 
0.506 0.041 


0.505 0.040 
0.501 0.018 


0.501 0.028 
0.500 0.014 


Ms 


0.825 0.280 
0.640 0.350 


0.565 0.190 
0.555 0.216 


0.532 0.142 
0.528 0.164 


(0.75,0.25) 


Vis 

i/ res 


0.774 0.121 
0.755 0.056 


0.755 0.052 
0.752 0.023 


0.751 0.036 
0.750 0.016 


Ms 
A res 


0.313 0.093 
0.300 0.146 


0.266 0.069 
0.268 0.094 


0.259 0.056 
0.260 0.072 


(0.95,5) 


i/ res 


0.969 0.131 
0.955 0.055 


0.953 0.058 
0.950 0.024 


0.952 0.042 
0.950 0.018 


Ms 

Ay-es 


11.251 3.492 
5.978 2.544 


5.836 2.104 
5.375 1.635 


5.397 1.645 
5.206 1.317 



Table 4.2 Average 95% confidence intervals for different values of v and A. 



(i/, A) 


Estimator 


n = 15 


n = 30 


n = 100 


ra = 500 


(0.1,1) 


^ls 

i/ res 


(0.060 , 0.158) 


(0.071 , 0.137) 


(0.063 , 0.142) 
(0.083 , 0.118) 


(0.084 , 0.117) 
(0.092 , 0.108) 


Ms 
\ res 
A res 


(0.214 , 53.963) 


(0.243 , 16.667) 


(-35.8067 , 58.864) 
(-0.362 , 2.885) 
(0.297 , 5.560) 


(0.108 , 2.228) 
(0.183 , 2.038) 
(0.464 , 2.652) 


(0.25,0.1) 


Vis 

l/ reS 


(0.151 , 0.389) 


(0.211 , 0.298) 


(0.162 , 0.359) 
(0.209 , 0.296) 


(0.211 , 0.292) 
(0.231 , 0.269) 


Ms 

X res 
A* 

res 


(0.015 , 2.183) 


(0.020 , 1.271) 


(0.034 , 0.201) 
(-0.031 , 0.291) 
(0.028 , 0.551) 


(0.052 , 0.159) 
(0.018 , 0.202) 
(0.044 , 0.265) 


(0.5,0.5) 


Vis 
V Tes 


(0.319 , 0.749) 


(0.366 , 0.664) 


(0.336 , 0.704) 
(0.424 , 0.586) 


(0.427 , 0.580) 
(0.464 , 0.536) 


Ms 

A res 
A* 


(0.099 , 9.154) 


(0.113 , 5.261) 


(-0.261 , 1.881) 
(-0.103 , 1.389) 
(0.164 , 2.542) 


(0.151 , 0.968) 
(0.121 , 0.978) 
(0.241 , 1.246) 


(0.75,0.25) 


Vis 

v res 


(0.527 , 1.057) 


(0.573 , 0.953) 


(0.534 , 1.019) 
(0.648 , 0.857) 


(0.648 , 0.855) 
(0.704 , 0.798) 


Ms 
X res 

A res 


(0.054 , 3.070) 


(0.067 , 2.206) 


(0.062 , 0.543) 
(-0.008 , 0.618) 
(0.087 , 1.055) 


(0.117 , 0.409) 
(0.076 , 0.453) 
(0.127 , 0.573) 


(0.95,5) 


Vis 
l/ res 


(0.719 , 1.221) 


(0.779 , 1.150) 


(0.700 , 1.219) 
(0.848 , 1.059 


(0.839 , 1.067) 
(0.904 , 0.999) 


Ms 
X res 
A* 

res 


(1.573 , 49.872) 


(1.602 , 30.456) 


(-2.5731 , 17.648) 
(0.285 , 10.951) 
(2.006 , 16.841) 


(0.963 , 10.334) 
(1.976 , 8.557) 
(2.714 , 10.067) 
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probabilities are closer to 95%, and has narrower in- 
tervals. Overall, the coverage probabilities and interval 
widths still provide good merits for our estimators even 
when the sample size is as small as n = 15. 

Collectively, Tables 4.1-4.3 strongly indicate that 
the proposed point and interval estimators performed 
well in our computational tests. We emphasize that 
the point estimates could also be regarded as reason- 
able starting values for better iterative estimation algo- 
rithms. 

5 Application 

We now apply our proposed methods to a real dataset. 
In particular, we estimate the parameters of the frac- 
tional Yule model using the branching times for plethod- 
ontid salamander dataset from Highton and Larson 
(1979) (see also Nee et al. 1994a; Nee 2001). The 25 
data points are the times measured from each node to 
the present of a phylogenetic tree, and can be down- 
loaded from the package laser of the R software. The 
summary statistics of the inter-branching times of the 
plethodontid dataset are given in Table 5.1 below. 

The point and the 95% confidence interval estimates 
are given in Table 5.2. The LS-based point estimate 
(0.749) of the fractional parameter v seemed to sug- 
gest that the plethodontid salamandar branching pro- 
cess is not a standard Yule process while the residual- 
based point estimate (1.119) appeared to suggest other- 
wise. Moreover, both the LS- and residual-based inter- 
val estimates of v indicated that v could be strictly less 
than one, which implies that a non-standard Yule pro- 
cess could model the plethodontid salamandar dataset 
with a confidence level of 95%. The residual-based point 
estimate (0.011) of A is more conservative than the 
bootstrap- and LS-based estimate (0.049). A similar 
observation can be gleaned from the 95% interval es- 
timates, i.e., the residual-based 95% interval estimate 
is narrower than the bootstrap- and LS-based interval 
estimates. 



We also tested the residuals for normality using the 
Shapiro- Wilk, Anderson-Darling, Cramer-von Mises, Lil- 
liefors, Pearson chi-square, and the Shapiro-Francia tests, 
which gave the p- values 0.811, 0.651, 0.619, 0.609, 0.849, 
and 0.461, correspondingly. Hence, these p- values in- 
dicated good fit of the fractional Yule process to the 
plethodontid salamandar data. 



6 Concluding remarks 

We have proposed closed-form expressions of the es- 
timators of the parameters v and A for the fractional 
linear birth or Yule, the fractional linear death, and 
the fractional sublinear death processes. The estimators 
were derived by taking advantage of the known struc- 
tural form of the logarithm of the random inter-event 
times and the well-studied least squares regression pro- 
cedure. The explicit formulas led to computationally 
simple and fast parameter estimation procedures. The 
inter-death time distributions and variances of the frac- 
tional linear and sublinear death processes were also 
obtained. These statistical properties were necessary 
for generating sample trajectories and for our estima- 
tion procedures to be applicable in these processes. It 
has also been shown that the proposed procedure can 
be easily extended to certain models that have differ- 
ent model parameterizations than the linear ones. The 
proposed methods were used to model a real physical 
process. Generally, the extensive computational tests 
showed favorable results for the proposed estimators. 

We cite some extensions which would be worth pur- 
suing in the future. For instance, improving the small- 
sample performance of the least squares-based estima- 
tors and developing other estimators using the likeli- 
hood approach or a re-sampling technique would be 
valuable pursuits. The application of these methods in 
practice, and the characterization of the appropriate 
functions m(9) and q(j) would also be of interest. 
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Table 4.3 Coverage probabilities and mean widths of 95% interval estimates for different values of v and A. 



(v, X) 


Estimator 


n = 
Coverage 


15 
Width 


n = 
Coverage 


30 

Width 


n — 
Coverage 


100 
Width 


n — 
Coverage 


500 
Width 




Vis 










0.949 


0.078 


0.952 


0.033 


(0.1,1) 




0.956 


0.098 


0.958 


0.066 


0.956 


0.035 


0.957 


0.016 


Ms 

A res 










0.893 
0.900 


94.671 
3.248 


0.915 
0.925 


2.120 
1.854 




A res 


0.919 


65.411 


0.928 


16.424 


0.946 


5.131 


0.946 


2.161 




V~ls 










0.937 


0.197 


0.958 


0.081 


(0.25,0.1) 




0.947 


0.238 


0.953 


0.162 


0.950 


0.087 


0.956 


0.038 


Ms 

X Tes 










0.940 
0.895 


0.167 
0.322 


0.949 
0.925 


0.107 
0.184 




A* 

res 


0.911 


2.168 


0.938 


1.251 


0.952 


0.521 


0.951 


0.219 




Vis 










0.946 


0.367 


0.933 


0.153 


(0.5,0.5) 




0.952 


0.430 


0.953 


0.298 


0.954 


0.161 


0.949 


0.072 


^ls 

X res 










0.931 
0.920 


2.141 
1.492 


0.923 
0.921 


0.817 
0.857 




A* 

res 


0.946 


9.055 


0.950 


5.148 


0.947 


2.379 


0.948 


1.004 




Vis 










0.931 


0.485 


0.942 


0.207 


(0.75,0.25) 




0.921 


0.529 


0.945 


0.379 


0.934 


0.209 


0.956 


0.094 


^ls 

X res 










0.941 
0.903 


0.481 
0.627 


0.950 
0.934 


0.291 
0.377 




A* 


0.916 


3.015 


0.931 


2.139 


0.951 


0.956 


0.952 


0.442 




Vis 










0.923 


0.519 


0.945 


0.228 


(0.95,5) 


Vres 


0.901 


0.502 


0.927 


0.317 


0.941 


0.210 


0.947 


0.095 


hs 
Ares 










0.899 
0.943 


20.221 
10.666 


0.939 
0.944 


9.371 
6.581 




A* 

res 


0.919 


48.299 


0.945 


28.854 


0.951 


14.748 


0.948 


7.356 



Table 5.1 Summary statistics for the plethodontid dataset. 



Minimum 


First quartile 


Median 


Mean 


Third quartile 


Maximum 


Standard deviation 


0.090 


0.607 


1.315 


3.981 


4.405 


22.120 


5.966 



Table 5.2 Point and 95% interval estimates for v and A of the plethodontid salamander data. 



Estimator 


Point estimate 


Interval estimate 


Vis 


0.749 


( 0.182 , 1.317 ) 


i^res 


1.119 


(0.955 , 1.283 ) 


A; s 


0.049 


( 0.008 , 0.089 ) 


A^es 


0.011 


( - 0.005 , 0.027 ) 


A* 

res 




( 0.003 , 0.051 ) 
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